Approximations for many-body Green's functions: insights from the fundamental 

equations 
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Several widely used methods for the calculation of band structures and photo emission spectra, 
such as the GW approximation, rely on Many-Body Perturbation Theory. They can be obtained 
by iterating a set of functional differential equations relating the one-particle Green's function to 
its functional derivative with respect to an external perturbing potential. In the present work we 
apply a linear response expansion in order to obtain insights in various approximations for Green's 
functions calculations. The expansion leads to an effective screening, while keeping the effects of the 
interaction to all orders. In order to study various aspects of the resulting equations we discretize 
them, and retain only one point in space, spin, and time for all variables. Within this one-point 
model we obtain an explicit solution for the Green's function, which allows us to explore the structure 
of the general family of solutions, and to determine the specific solution that corresponds to the 
physical one. Moreover we analyze the performances of established approaches like GW over the 
whole range of interaction strength, and we explore alternative approximations. Finally we link 
certain approximations for the exact solution to the corresponding manipulations for the differential 
equation which produce them. This link is crucial in view of a generalization of our findings to the 
real (multidimensional functional) case where only the differential equation is known. 

PACS numbers: 71.10.-w 



I. INTRODUCTION 

The one-particle Green's function (GF) 1-3 is a pow- 
erful quantity since it contains a wealth of information 
about a physical system, such as the expectation value 
of any single-particle operator over the ground state, the 
ground-state total energy, and the spectral function. In 
order to access this quantity one can start from its equa- 
tion of motion: 4-6 



In (3) Go determines the appropriate initial condition 
in time; note that the solutions of (1) and (2) are not 
unique. Moreover, in order to calculate G, the knowledge 
of G 2 is required (which in turns requires the knowledge 
of G3, and so on) 4,6 . In order to obtain a closed expres- 
sion one can generalize G(l,2) to G(l,2; [<p]), where an 
external fictitious time-dependent potential ip is applied 
to the system. This allows one to express G2 as 7 



' 9 U ^ 



G 2 (3,4;2,4+;M) 



G(l,2) 



G(3,2;M)G(4,4+;M) 
SG(3,2;[tp]) 
M4) ' 



(4) 



+ 1 / (i3u(l+,3)G 2 (l,3;2,3 + ) = <S(l,2), 



G) 



where h(r\) is the one-electron part of the many-body 
Hamiltonian, G 2 (l,3;2,3 + ) is the two-body Green's 
function, and i;(l + ,3) is the Coulomb potential. The 
space, spin and time variables are all combined in (1) = 
(n, oi, ii), and (1 + ) = (n,<7i,ii) with tf = t\ + 6 
(<5->0+). 

Eq. (1) can be manipulated in order to get a more 
practical expression, by introducing the non-interacting 
Green's function Go with 

- h(n)\G (l,2) = 6(1,2), (2) 

which reinserted in Eq. (1), gives 
G(l,2) = G (1,2) 

-i J d3d4G (l,3)w(3 + ,4)G 2 (3,4;2,4+). (3) 



Note that in (4) all Green's functions are generalized to 
non-equilibrium since they depend on the perturbing po- 
tential. The equilibrium G and G 2 in (3) are then ob- 
tained by taking <p = 0. Inserting (4) into (3) yields a set 
of functional differential equations 4 for the unknown G 

G(1,2;M) =G (1,2) 

d3 Go(l,3)V- ff (3;[^])G(3,2;M) 



d3 G (1,3V(3)G(3,2;M) 

£ G(3,2; [<p]) 
M4) : 



+ i / d4d3G (l,3)t;(3+,4)- 



(5) 



where Vh(3) 



dAv{3, 4)G(4, 4+; [ip]) is the Hartree 



potential. Since the Hartree potential contains the 
Green's function, this term makes the equations nonlin- 
ear. We are interested in the solution of Eq. (5), for 
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tp = 0. Its calculation would hence require the solution of 
a set of coupled, non-linear, first-order differential equa- 
tions, which is clearly a non trivial task. Moreover one 
would need a new initial condition to completely define 
the desired solution of this differential equation, since 
SG 

the derivative — — has been introduced. Therefore usu- 

Oifi 

ally another route is taken: one includes the functional 
derivative in (5) in the definition of a self-energy 4 



£(1,3) = i J d4d2v{l + A) 
x G- 1 (2,3), 



5G(1,2;M) 



M4) 



<p=0 



(0) 



which, inserted into Eq. (5) for tp = 0, gives: 

G(l, 2) = G (l, 2) + / d3 G (l, 3)^(3)G(3, 2) 



d4d3 G (1,3)£(3,4)G(4,2) 



(7) 



This is the Dyson equation for G, where £ contains all 
the many-body effects (beyond the Hartrcc contribu- 
tion) present in the system. Of course |^ and therefore 
£ are still not known and, in practice, £ has to be 
approximated. A good starting point is obtained by 
reformulating the problem in terms of a coupled set of 
equations containing the one-particle Green's function, 
the polarizability P. the self-energy £, the screened 
Coulomb interaction W, and the vertex T. These 
equations are most often solved within the so called GW 
approximation (GWA) 8 , where the vertex T is set to 
unity, resulting in £ w iGW. Over the last two decades 
the GW method has become the tool of choice for 
calculations of quasi-particles (QP) band structures 9 ' 10 
of many materials and direct and inverse photo emission 
spectra (see e.g. Rcf. n ~ 14 ) improving substantially 
over the results provided by static mean-field electronic 
structure methods. 

However the GW A suffers from some fundamental 
shortcomings (see e.g. Refs. 15-18 ) and with £ being 
of first order in W, is not expected to describe strong 
correlation. Higher orders in W could be added by 
iterating the equations, but this is technically difficult, 
and there is no guarantee that results will quickly 
improve. It is therefore necessary to find guidelines. 
In the present work we go back to Eq. (5). Our 
aim is first, to obtain new insights about standard 
approximations by relating them more directly to the 
original equations. Second, we want to use Eq. (5) to 
explore alternative approximations. Finally, it might be 
interesting to concentrate directly on the set of coupled, 
non linear, first order functional differential equations 
for G, Eq (5), although it has been acknowledged that 
no "practical technique for solving such functional 
differential equation exactly" 4 is available. However, one 
may still hope that with new algorithms and the increase 
in computer power, numerical solutions might become 
accessible. The present work is hence also meant to 



explore strategics for, and possible problems of, such a 
route. In the following we resort to two approximations. 
First we linearize the set of equations by expanding Vh 
in terms of tp. Second, we discretize Eq. (5) and consider 
in a first instance only one point for each space, spin, 
and time variable: we will call this latter approximation 
the "1-point modeV\ as opposed to the full functional 
problem. The strategy underlying this procedure is 
the following: for the 1-point model, we can derive the 
exact explicit solution of the now algebraic differential 
equation, and solve the initial value problem. One can 
hence explore approximations to the full solution, which 
yields valuable insights in the performance of current 
approaches and suggestions for alternative ones. By 
determining which manipulations of the differential 
equation (DE) produce such approximate solutions, one 
obtains suggestions for analogous manipulations on the 
differential equation for the full functional problem, 
which opens the way to translate our model findings to 
realistic calculations. 

The manuscript is structured as follows. In Section II we 
present the linearized differential equation which can be 
solved exactly within the 1-point framework. We discuss 
in particular the initial value problem, and how it can 
be overcome. In Section III we examine, in the 1-point 
framework, various common approximations to the solu- 
tion of the DE: the iteration of Eq. (5) and approxima- 
tions based on a Dyson equation, in particular different 
GW flavours. In Section IV we explore other routes to 
manipulate the initial DE and obtain approximate solu- 
tions. We finally give our conclusions and perspectives 
on future work in Section V. 



II. THE SCREENED EQUATION IN A 1-POINT 
FRAMEWORK 



Our first goal is to simplify the equations such that the 
main physics is retained, but manipulations become more 
straightforward. To this end we linearize the differential 
equation with an expansion of the Hartree potential to 
first order in the external potential tp, 



V H (3; [tp]) « -i I d4«(3+,4)G(4,4+; [tp]) 
SG(4A+;[tp}) 



d4d5v(3 + A)- 



M5) 



¥J=0 



ip=0 

tp(5)+o(tp 2 ).(8) 



Eq. (5) hence becomes 
G(l, 2; [?]) = G° H (1,2) + I d3G° H (l, 3)$>(3)G(3, 2; [<p]) 



, d3d5G°,(l,3)IF(3 + ,5)^|M ;(9) 

where G a H is a Hartree Green's function containing the 
Hartree potential at vanishing tp, tp = e~ 1 tp is the renor- 
malized external potential, and W = e~ 1 v is the screened 
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Coulomb potential with e the dielectric function at tp = 
19 . 

Concerning Eq. (9) three important remarks should be 
made: first, through the linearization the screened in- 
teraction W becomes the central quantity of the equa- 
tion. This is justified by the physics of extended systems, 
where screening and plasmons are key concepts. 
Second, in principle W is the exact screened interaction, 
which of course is not known. One can however adopt 
two strategies: either W is considered to be an exter- 
nally given quantity, obtained within a good approxi- 
mation, e.g. from a time dependent density functional 
theory (TDDFT) calculation 20 ; or one could also recal- 
culate W from G[0\ (see in the next section). In this 
work we will adopt the first strategy, which is illustrated 
in Fig. 1. Such a philosophy is rigorously justified. In 
particular, in the framework of the theory of functionals 
it is possible to pass from the Luttinger-Ward functional 
(given as functional of G, though indeed one should add 
the bare Coulomb interaction v as argument) to the so- 
called ^-functional, where v is replaced by W 21 . This 
explains for example why self-consistency in G only (and 
not in W) is sufficient to have a conserving GW approxi- 
mation. Moreover, in practice this is the most current 
way of proceeding, corresponding e.g. in a GW cal- 
culation to the "best G, best W" approach: while the 
non-interacting G is taken e.g. to be the Kohn-Sham 
Green's function, W is calculated as accurately as possi- 
ble, e.g. in the adiabatic local density approximation to 
TDDFT. Third by approximating the functional deriva- 

5G SG^ 1 
tive — = —G G ps GG (which supposes the self- 

oip dip 

energy to be independent of (p) one obtains the Dyson 
equation for the one-particle Green's function in the GW 
approximation to the self-energy. The proof is given in 
App. A. This result shows that, even though the lin- 
earization procedure is an approximation, Eq. (9) is still a 
promising starting point to analyze the different flavours 
of the GW approximations and to go beyond. 
After linearizing, the next step consists in discretizing 
Eq. (9) and then in considering only one value for the 
space, spin, and time variables respectively (or cquiv- 
alently concerning space and spin, in considering all 
Green's functions to be diagonal in a given basis): this 
is the 1-point model employed throughout the whole 
manuscript. The 1-point framework has already been 
used by other authors: in Refs 22 ' 23 Hedin's equations 
are combined in one single algebraic differential equation 
which is solved as a series expansion. This allows the 
authors to enumerate the diagrams for a certain order of 
expansion. Several expansion parameters are examined, 
such as vG'jj, with v the bare Coulomb potential and 
Gh the Hartree Green's function, vG 2 , with G the exact 
Green's function, WG 2 , with W the screened Coulomb 
potential, etc., which shows how at various orders of ex- 
pansion the number of diagrams decreases by increasing 
the degree of renormalization. This is also the spirit be- 
hind the linearized equation (9) , in which the natural ex- 



pansion parameter would be WG 2 H: where W is treated 
as an externally given interaction. The advantage of us- 
ing the 1-point framework is that the equations become 
algebraic and thus the enumeration of diagrams is fa- 
cilitated. In Ref. 24 a similar strategy as in Refs 22,23 is 
used to enumerate diagrams, focusing in particular on the 
asymptotic behavior of the counting numbers. Moreover 
Hedin's equations are transformed into a single first oder 
differential equation for the GF as a function of an inter- 
action parameter, and an implicit solution is obtained. 
In order to fix the particular solution of this differential 
equation the initial condition Gr v =o) = Go is used. In- 
stead here we concentrate on (5), or better its linearized 
form (9), which is another differential equation for G, as 
a functional of an external potential. This choice allows 
us to (i) emphasize the essential physics contained in the 
screened Coulomb interaction W, (ii) discuss various as- 
pects of the many-body problem in a clear and simple 
way, (iii) obtain an exact solution of the approximate 
equation that can be used as a benchmark. Moreover we 
believe that the 1-point version of Eq. (9) can be a natural 
starting point for a generalization to the full functional 
problem. While the equations are easier to manipulate, 
physical information is of course lost in the 1-point frame- 
work. In particular no poles (addition/removal energies) 
of the GF appear. However, the various aspects that will 
be explored in the following are intrinsically related to 
the structure of the equations, and hence exportable also 
to the full functional problem, in the same spirit as in 
Refs 22 " 24 . 



A. The 1-point differential equation 

In the 1-point model Eq. (5) reduces to an algebraic, 
non-linear, first order differential equation 

y{z) = Vo + vy y 2 (z) + yozy(z) - vy ^ ^ * (10) 

where (p -> z, G(l,2; [<p]) -> y(z), and G (l,2) -> y . 

Moreover ii>(3 + ,4) — > — v: this change of prefactor com- 
pensates for the time- or frequency integrations that have 
been dropped in the 1-point model and corresponds to a 
standard procedure in this context 22,24 . We can now lin- 
earize Eq. (10) in the same way as we did starting with 
Eq. (5) and obtaining Eq. (9). This yields 

y u {x) = y H + yHxyu(x) - uy H ———. (11) 

Hence with respect to Eq. (9), Cp — > x, G° H (1,2) — > y$ , 

and iW(3 ,5) — > —u: the subscript u in y u highlights 
its u dependence. In the following, for simplicity of no- 
tation, we denote y° H by yo unless stated differently. In 
Appendix B we sketch the main steps to solve Eq. (11), 
based on the general ansatz y u (x) = A(x) ■ X(x). With 
the choice 

A(x) =e[ fe "^] (12) 
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one obtains the equation 

dl(x) 1 



d.r 



—e 

u 



(13) 



and the general solution y u (x) reads 
Vu{x) 



x { erf 



I 1 \ 1 

(v )\ 7T 

Ua V 2u 



C(y ,u)}, (14) 



where G (yo j w ) is to be set by an initial condition. In 
the limit x — > 0, which is the equilibrium solution we are 
looking for, Eq. (14) becomes 



g*»y 



x { erf 



1 



C{yo,u) 



(15) 



2uyl\ 

Note that a similar ansatz can also be used for the 
full functional problem, namely G(l,2) = J d3A(l,3) • 

1(3, 2), in order to get a set of differential equations that 
are less complicated to manipulate than the original one. 



B. The initial value problem 

In general in order to set C(yo,u), y u (x) has to be 
known for a given potential xp (i.e. y u (xp) = j/f ). How- 
ever it is far from obvious to formulate such a condition 
in the realistic full functional case; this would indeed re- 
quire the knowledge of the full interacting G for some 
given potential ip. Therefore the question is whether one 
can reformulate the condition in a simpler way in order 
to set G. 

To answer this question we expand the exact solution for 
small values of u, obtaining: 



°(l + C(u,y ) 



+ 2/o 



uyl 



iu 2 yl - lbu z yl + o(u 4 )Y (16) 



When u — > the one-body Green's function G has to 
reduce to the non interacting Go, in our framework this 



translates into: y l 
Eq. (16) gives 



2/o- Imposing this condition on 



2 a 

which is satisfied if 



^\(l + C(u,y j) =0, u->0, (17) 



C{u,y ) 



-1. 



u -> 0. 



(18) 



This result for C holds also for u ^ 0. Indeed it guaran- 
tees a non-divergent result for any non-vanishing poten- 
tial x in (14). Moreover it reproduces the perturbativc 



result which is obtained by iterating Eq. (11); for exam 
pie the sixth iteration yields 

uyl + Zu 2 yl - lbu z yl. 



yo 



(19) 



This is precisely the same series as the one appearing 
in Eq. (16) when C(u,yo) is set to —1. 



III. ANALYSIS OF COMMON METHODS TO 
CALCULATE THE ONE-BODY G 

In the following we will analyze various established ap- 
proximations for the calculation of the one particle G, 
using the knowledge of the exact solution. 



Iteration of the DE 



Let us first iterate Eq. (11) starting from y^\x 



according to: 



= yo + y { u ] xyo - uito^H (20) 

ax 



For x = the first two orders in u read 



y { u ] = yo - uyl, 

y { u ] =Vo- uyl 



3u 2 y 5 , 



(21) 
(22) 



and Eq. (19) for the third order. Results as a function 
of u are depicted in Fig. 2 together with the exact solu- 
tion. Two observations can be made: i) very few terms 
are needed to obtain a good approximation to the exact 
solution in the small u regime; ii) for a given u — u n , the 
expansion diverges starting from an order n. The larger 
is u n , the smaller is n, which limits the precision that can 
be obtained. As previously mentioned, the iteration coin- 
cides with the expansion for small u of the exact solution. 
Since the small u expansion is de facto the asymptotic 
expansion of the error function times an exponential (as 
can be seen in (16)) the divergent behaviour of the it- 
eration in (20) is not surprising. Divergences of higher 
orders have been found in perturbation expansions for 
realistic systems, e.g. for orders higher than 3 in the 
M0eller-Plesset scheme 25,26 . 



B. Self-energy based approximations 

In this section the introduction of a self-energy E will 
be discussed along with its most common approxima- 
tions. 

The Dyson-like form for Eq. (11), which is the equivalent 
of Eq. (7), reads: 

y u (x) = t/o + yoxyu(x) + y ^ u [y u (x)] y u (x) (23) 

where a self-energy kernel has been defined as 

dy u (x) 1 



S„ [y u {x)} = -u- 



dx y u (x)' 



(24) 
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l ^4^- = -yt(x) dy \^ ■< 



dx 

r„ [y u (x)] = - 

energy reads 



dx 



dx 

for the vertex function, the self- 



S u [yu(x)] = -uy u (x)T u [y u (x)] 



(25) 



which is the equivalent of £ = iGWT s . The Bethe- 
Salpeter equation for the vertex function V is then de- 
rived from (23) 



dy u 1 (x) 
dx 



■1 



dT, u [y u {x)\ 
dx 

dT, u [y u (x)} dy u (x) 



dy u (x) 



dx 



(26) 



We call this procedure iterative self-consistent scheme, in 
contrast with the direct self-consistent scheme where one 
solves directly the Dyson equation (23), for x = 0, with 
= —uy u . In this latter case one gets a second-order 
equation with two solutions 



±Vl + 4uyg-l 
2uy 



(32) 



Note that for the full functional problem one would find 
even more solutions, since a second-order equation has 
to be solved for each matrix element of G. 
In order to choose the physical solution we Taylor expand 
the square root around u = 0, which leads to 



from which for x — > 



T u (y u ) = l + ^yAr u (y u )yl (27) 
dy u 

where y u = y u (x -> 0). For x = Eqs. (23), (25), 
and (27) correspond to a subset of the so-called Hedin's 
equations 8 , obtained by fixing W. A pictorial represen- 
tation of this subset for a given W is given in Fig. 1. 
In the following we will approximate the equations and 
the results will be compared to the exact solution of the 
differential equation, in order to obtain greater insight 
about these self-energy based techniques. From now on 
all quantities will hence be understood to be taken at 
x = 0. 



1. GoWo and self-consistency 

Let us first look at different flavours of the GW 
approximation 8 . Setting T u (y u ) to unity, it follows that 

S«(y«) = —uy u . Within the initial guess y^ = yo, one 
obtains a so-called GqWo self-energy E u = —uy - 27 This 
is then employed in the Dyson equation (23) in order to 

get an improved y« . To go beyond this first approxima- 
tion one can iterate further within the GW approxima- 
tion, i.e. keeping r„ = 1. This corresponds to an iter- 
ation towards a GWo result, since G is iterated towards 
self-consistency but u, which represents the screened in- 
teraction, is kept fixed. We report here the expressions 
obtained for Go Wo, i.e. the first solution of the Dyson 
equation, and for three successive loops 



yi 1] 



!J 



G Wo _ 



!Jo 



yi 2) 



yl 3) 



yi 4) 



i 



= yo 



yo 



I/O 



uyl 



1 + 2uyl ' 

1 + 2uyg 
1 + Zuyl + u 2 y\ ' 
1 + iuyl 



2 4 



1 + 4uyg + 3u 2 y% ' 



(28) 
(29) 
(30) 
(31) 



yu~±[yo + 



l 



1 



2uy ) 2uy 



(33) 



Since for u = 
solution is y u 



one has to obtain y u 
= y/1 + Any 2 - 1 
2uy 



yo, the physical 
In Fig. 3 we can appre- 



ciate how well these GW-based methods are performing 
against the exact solution in a wide u range. 
Interestingly odd iterations quickly converge to the phys- 
ical solution of the direct sc-GW-o, while even iterations 
do also converge but at a slower pace: it can be shown 
that for u — > oo their limit forms the sequence of rational 



which ultimately converges 



numbers I—;—;—;—; — - 

\2' 3' 4' 5' 6 

to 0. All the odd iterations have instead the exact large 
u limit (namely y u = when u — > oo). One might use 
this property to improve the convergence of the series. 
An important question is now: does the result of the self- 
consistent procedure depend on the starting point of the 
iteration? Here we have naturally chosen y^ = y , but 
one might fear that this choice is simply lucky. Let us 
therefore look at the general iterative scheme which is 
obtained by solving the Dyson equation (23) for x = 



y u = 



Do 



1 + youy u 



(34) 



By starting the iteration with a guess for y u on the right 
side one obtains 



y { u +1) 



yo 



i , («)' 

1 + youy u 



(35) 



For y u ^ = y s one has e. g. after the third iteration 



yi 3) 



yo 



u Vo 



(36) 



l+y uy 3 



6 



This contains nothing else but the continued fraction rep- 
resentation for the square root 



ViTz~ = 1 



z/2 



z/A 



(37) 



z/A 



z/A 



1 + 



z/A 



corresponding to the physical solution y u 



y/T+Z-1 

2uy 



where z = Auy$. It converges for all values of the termi- 
nator y s . Therefore, this iteration will always converge 
to the physical solution. Does this mean that there is no 
risk of running into the unphysical solution? The answer 
is that it depends on the iterative scheme that is used, 
and not on the starting point. Look at the following way 

to re- write the Dyson equation (23): —uy u = (in 

2/o Vu 

other words, S = G^ 1 — G -1 ). If we iterate this equation 

by starting with some y^ = y s on the right-hand side 
we get 



i 

uyo 



1 



hence 



2uy y 



a y 



2uyl 



(n)- 



(38) 



(39) 



1+- 



!+ V 



which, with Eq. (37), is just the continued frac- 
tion represen tation for the unphysical solution y u = 
(— yl + Auvq — l)/2uyo- In a way, this is good news: 
usually the iterative scheme adopted in the context of 
GW calculations is rather the first, safe one. Indeed it 
has been found empirically that such a scheme leads to 
self-consistent results independent of the starting point 
and in reasonably good agreement with experiments (see 
e.g. 28-30 ). However, when one goes beyond GW, higher 
order equations appear as we will see in the following. 
There arc hence more and more solutions, and more and 
more ways to iterate the equations. In other words, there 
will be an increased danger to run into a wrong solution. 
One should keep this in mind when trying to add vertex 
corrections beyond GW. 



2. Vertex corrections - First order F 

We will now analyze the effects of a first order vertex 
correction which is obtained employing S u = —uy u in 
Eq. (27) 8 ' 31 . Solving for T u gives 



1 



(40) 



Employing this vertex the self-energy (25) becomes 

1 



^\yu) = - 



uy u 



1 



uyl 



(41) 



Now two routes can be taken and either a GoWoT^^yo) 
or a self-consistent GWqT^ (y u ) calculation can be car- 
ried out. The first of the two is once more based on 
the initial guess for the Green's function y^ = yo, and 
consequently the vertex and the self-energy in ( 40-41) 

read respectively T^(y ) = — ^ and T,£\y Q ) = 

1 + uy£ 

Solving the Dyson equation with the 



-uyo 



1 



1 



above ingredients yields: 



y 



GoWaT 



yo (l + uyf) 
l + 2uy 2 



(42) 



Instead, solving the Dyson equation in a self-consistent 
fashion, with the expressions (40-41) yields: 



,GW r 




1 



27u 3 



1 

Au 2 



27u 3 



5>- 



As it can be noticed from the result a cubic equation 
for the unknown y u had to be solved within this more 
sophisticated approach. Again the limit of vanishing in- 
teraction has been used to pick the physical solution. 
In Fig. 4 we can directly compare the two types of ver- 
tex corrections. For small u values their performance 
is similar, however, in a wider u range (see inset), the 
GqWoT^ scheme diverges from the exact solution and 
has the wrong asymptotic limit u — > oo: it hence behaves 
as the first iteration of the sc-GWo approach, which also 
exhibits the wrong large u limit. Fig. 4 also shows how 
the GWqT^ scheme, for small u values, slightly improves 
over the sc-GWo- However, given the augmented com- 
plexity already at this first order of the correction (one 
could very well iterate further the equations for T and £ 
and get higher order corrections), the benefits of employ- 
ing vertex corrections are not obvious. Also note that 
interestingly, on the scale from u — to u — > oo, the 
closest curve to the exact one is the sc- GWq one. 



IV. EXPLORING OTHER APPROXIMATIONS 
FOR G 



In this section we will explore alternative approxima- 
tions to the exact solution of the 1-point DE and the cor- 
responding manipulations of the initial differential equa- 
tion producing them. Here we will report in particular 
approximations that might be eventually transposed to 
the full functional framework. 
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A. Continued fraction approximation 

A well known approximation for the error function is 
its continued fraction representation 32 . The exact ex- 
pression for y u (Eq. (15)) transforms into 

1 



2u 



1/2 



(44) 



3/2 



\/ 2u Va 



I/O 



u Vo 



(45) 



2uyf ) 



1 



We will now show how one can obtain Eq. (45) starting 
simply from the initial DE in Eq. (11), equivalent to (9), 
without any information about its exact solution. Begin- 
ning with Eq. (11) and taking successively higher order 
derivatives of the equation, one obtains : 

dy u (x) , . dy u (x) d 2 y u (x) 

--yoy u (x) + Vox— uy (46) 



dx 
d 2 y u (x) 

dx 2 
d 3 y u {x) 

dx 3 



= 2w 



: 3y ( 



dx 
d 2 y u {x) 



dy u (x) 
dx 

d 2 y u {x) d 3 y u (x) 



Vox- 



dx 2 



dx 2 



-yox- 



dx 3 



dx 2 
d 3 y u {x) 
dx 3 

d i y u {x 
dx 4 



uyo- 



uyo- 



(47) 
48) 



and so on. Neglecting derivatives e.g. from the 4 th order 

on and then setting x = 0, this truncation allows us to 

solve all the above equations, beginning with Eq. (48) 

i ii- ■ d 3 y u (x) 

(now an algebraic equation in the unknown — — 5 — by 

dx 3 

d 2 y u (x) . 
keeping — — — — as parameter); subsequently we insert 

2 

the result in (47) and solve for ^ ^ £ Qr 



and ultimately Eq. (11) getting 



yo 



dx 2 



dx 



(49) 



1 + 3uy 2 



1 



which is precisely the result obtained by approximating 
the exact solution with a continued fraction expression 
for the error function ( Eq. (45) ). We will name this ma- 
nipulation limited order differential equation. In Fig. 5 we 
compare the different orders of this approximation to the 
exact expression for y u . The approximation gets rapidly 
closer and closer to the exact solution by including higher 



derivatives. However, also for this continued fraction, 
odd and even orders, converge towards the exact result 
with a different speed. In analogy with the continued 
fraction of Eq. (37), even iterations have the correct large 
u limit, while the odd ones don't, although they do even- 
tually approach it for a very large number of steps. We 
notice that the above continued fraction converges slower 
than the one arising from the sc-GWo; however, the for- 
mer will eventually converge towards the exact solution, 
whereas the latter only to the self-consistent GWq result. 
It is therefore interesting to note that such a procedure 
can in principle be used also in the full functional frame- 
work (see related manipulations e.g. in 6,33 ), where the 
functional differential equation can be differentiated to 
an arbitrary order and the corresponding approximated 
G obtained. For example, differentiating Eq. (9) with 
respect to the external potential Cp one gets 



<*?(!, 2; [£]) 
M6) 



rf3G H (l,3)^||G(3,2;[^) 
6G(3,2;[0\) 



d3G° H (l, 3)v>(3) 



+ i / d3d5M / (3 + ,5)G? f (l,3) 



^G(3,2;[^) 



(50) 



S 2 G 



Truncating the highest order derivative -r^r and solving 

0(f z 

for tp = (which means also (p = 0) gives: 



SG(1,2;[0\) 



= G° r (l,5)G(5,2) 



(51) 



^(5) 

which reinserted in Eq. (9) yields: 
G(1,2) = G°,(1,2) 

+ i J d3d5G° H {l,3)W(3 + ,5) 

x G^(3,5)G(5,2). (52) 

Like in the 1-point model, this first step simply provides 
the one-particle GF in the G^Wq approximation to the 
self-energy. One can go further: differentiating Eq. (50) 
with respect to (p and neglecting the third order deriva- 
S^G 
S>p 3 



five -ttt yields: 



G(l, 2) = G^(l, 2) - il d5d3d8d9G° H (l, 3)W(3+ 5) 
x m _1 (3 5 5 ; 9 = 8 )Gff(9, 8)G(8, 2) (53) 

with 

m(16;57) : = -5(15)5(76) 

+ i J d3G° H (l,3)W(3 + ,5)6(7,6) 
x[G° H (3,6) + G%(3,5)], (54) 
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which is a four-point quantity of a similar complexity as 
the Bethe-Salpeter equation 5 . Indeed in the full func- 
tional problem the equations become quite involved since 
terms like uy 2 correspond to large matrices. However 
the approach doesn't require self-consistency. This might 
turn out to be a significant advantage, compared to ver- 
tex corrections to E, as we have discussed in the previ- 
ous subsection concerning self-consistency. More details 
about the derivation are given in App. C. 



one obtains for the different orders of the full solution 



B. Large u expansions 



Perturbation theory usually deals with weak interac- 
tions, hence the small u limit. However, it is also very in- 
teresting to examine the large u limit for several reasons: 
i) this is the regime of strong correlation, where current 
approximations exhibit failures; ii) the large u expansion 
of the exact solution gives a convergent series (being a 
product of two convergent Taylor expansions, one for the 
exponential and the other one for the error function) and 
one can, for instance, obtain a better approximation to 
the exact solution by adding higher order terms (which 
instead does not improve the result for the small u ex- 
pansion of the solution) ; iii) excellent approximations for 
the exact solution are Pade approximants 34 , which have 
to be constructed using both the small and large u limit. 
In this subsection we will present two possible routes to 
approach this limit: the first is a straightforward large u 
expansion of the exact solution for y u , while the second 
combines the latter with the large u expansion for the 
Dyson equation. 



y ( u 1/2) = 



yi~ 3/2) 



yi~ 2) - 



yt 5/2) = 



+ 



n 
2u 
1 

uyo 
1 

uyo 
1 

uyo 
1 

uyo 
1 

]u 2 y$ 



1 



yi~ 3) = 



i 

— h 

uyo 
1 

10u 3 j/q 



2uy 2 
1 

2uy 2 
1 

2uy 2 

IT 

2u 



2uy 2 



TC 

2u 

7T 

2~77 

7T 

2u 
2u 

7T 

2~77 



TT 

277 
l 



1 

§u 2 yl 



1 



1 



§u 2 yl 8u 2 y 4 



(57) 
(58) 
(59) 
(60) 

(61) 

7T 

277 

(62) 



Fig. 6 shows how these different expansions perform ver- 
sus the exact result. Overall their behaviour is very good 
for large u and few orders are sufficient to get a good ap- 
proximation over a wide u range (which is our ultimate 
goal), however for u = they all diverge. 



2. Large u expansion for y u and for the Dyson equation 

When u gets larger, also S„ increases. This im- 
plies that, using the Dyson equation for the one-particle 
Green's function y u = (yj^ 1 — S u ) one could expand 
y u as 



Hence to lowest order y 



(63) 



-S: 1 or 



1. Straightforward large u expansion for y u 



By expanding both the exponential prefactor and the 
error function appearing in Eq. (15): 



Vu 



(64) 



This simple relation allows us to use the large u expansion 
of the exact solution for y u to approximate E u for large 
u; we can then use this approximate £„ in the Dyson 
equation to recalculate y u . For example, using the lowest 
order of the large u expansion of the exact y u one gets 
the following self-energy: 



erf 



2uy$ 



2uy 2 8u 2 y 4 



2 r 



1 



1 



1 



40u 2 y% V 2uy 2 



(55) 



1 



2uy^ 6uy$ V 2uy£ 



(56) 



(65) 



which reinserted in the Dyson equation y u = (y 1 



gives: 



Vu 



yo 



(66) 
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In Fig. 7 the performance of this approximation for y u is 
plotted against two orders of the straightforward large u 
expansion for the Green's function, GqWq and the exact 
solution. The " large £" approach shows an overall good 
agreement (generally better than Go Wo) with the exact 
solution and it has the desirable property of being exact 
in the small and large u limits, mending the divergence 
of all orders of the straightforward expansion for y u . At 
higher orders of the approximation this property remains 
true, although undesired poles appear. In conclusion the 
methodology is promising and worthwhile to be explored 
further. The main difficulty is that in the framework of 
a large u expansion, without knowing the exact solution, 
one would not straightforwardly know how to set the con- 
stant G, i.e. how to pick the physical solution: this issue 
requires further analysis. 



C. Self-consistent calculations of the Hartree 
Green's function and of the screened interaction 

In the above discussions we have treated the Hartree 
Green's function and the screened interaction as exter- 
nally given quantities. This is justified by the fact that 
realistic calculations are most often following such a prag- 
matic ansatz. In principle these quantities should be part 
of Hedin's self-consistent cycle. A fully self-consistent 
treatment, in the full functional framework, is today out 
of reach. In the 1-point model, however, it is possible to 
go beyond this limitation and indeed, the implicit solu- 
tion of Hedin's equation that has been achieved in the 
work of 24 contains all quantities calculated on the same 
footing. Also in the linearized version, that is employed 
in the present work, one can obtain the Hartree Green's 
function and the screened potential consistently from the 
equations, as we will discuss in the following. Let us first 
turn to the Hartree Green's function i/jj . In terms of the 
truly non-interacting Green's function yo it reads 



/In 



youy u 



(67) 



in other words, it depends (through the density) on the 
solution y u at vanishing external potential. In a self- 
consistent scheme this y® H should then replace yo in the 
solution Eq. (15), which leads to an implicit equation for 
y u . For a self-consistent treatment of the screened inter- 
action we can use of the fact that the 1-point differential 

equation can be solved for — ^ , and insert the result into 

ax 

the expression for the screened interaction u in terms of 

the bare v, which reads u = v + v—^v. Two routes can 

ax 

be taken. The first one is based on the linearized equa- 
tion (11) where the interaction is already screened from 
the very beginning. This leads to a quadratic equation 
for u, with two solutions 



(68) 



The physical solution is the one of the positive square 
root, since it approaches the bare v in the limit of van- 
ishing interaction, hence vanishing screening. The second 

route consists in calculating — ^ from the initial equa- 

dx 

tion (10), where the bare yo and the interaction v appear. 
This yields 



Vu 

u = v I 2 



vyl 



(69) 




In both cases, the solution for u should be used into 
Eq. (15), which again makes the expression for the GF 
implicit. One may argue about which of the two ways 
to calculate u sclf-consistcntly is more adequate. In a 
realistic calculation one would probably use the former 
approach, in an iterative way: after calculating the GF as 
a functional of the external potential for a given initial in- 
teraction in the linearized DE, one would recalculate the 
W from the functional derivative, and so on. Whatever 
choice, however, does not influence the main conclusions 
that can be drawn from the above considerations. Specif- 
ically: i) a self-consistent calculation leads to an implicit 
solution (like in the work of 24 ) which however would not 
be identical to theirs because of our linearization proce- 
dure; ii) the behaviour for the small interaction limit is 
unchanged by the self-consistent treatment, as one can 
verify from equations (67), (68) and (69); this means in 
particular that the constant G is chosen in the same way 
as before, iii) Finally also the discussion about the limit 
of large interaction is unchanged: by making the ansatz 

that to lowest order y u oc — = one finds consistency. 

Altogether, this shows that the linearization of the equa- 
tions docs not imply necessarily that one has to treat 
the Hartree Green's function and the screened interac- 
tion as externally given quantities. It also shows that a 
more refined, self-consistent treatment does not change 
the overall behaviour of the solution. 



V. CONCLUSIONS AND OUTLOOK 

In this paper we explore several aspects of the set of 
first order nonlinear coupled differential equations which 
are conventionally solved perturbatively in order to cal- 
culate the one-particle Green's function. After the lin- 
earization of the Hartree potential with respect to the 
external one, we employ a 1-point model where the set 
of -now linear- differential equations reduces to a 1 st or- 
der algebraic differential equation, that can be solved ex- 
actly. This provides insights into the structure of the gen- 
eral family of solutions, and on how to determine which 
amongst them corresponds to the physical one. Within 
the model we study the performance of established ap- 
proaches over the whole range of interaction strength: we 
find that iterations towards self-consistency in the GW 
scheme sensibly improve on the one-shot (Go Wo) calcu- 
lation and that including first order vertex corrections 
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improves the self-consistent GWo results only slightly 
and only for small u. We also find that in case of self- 
consistent GWo two solutions are possible, of which only 
one is physical. We show that the standard iterative 
scheme will always converge to the physical solution, al- 
though other schemes may yield different results. This 
is an important finding: when going beyond GW both 
the number of possible solutions for the Green's function 
and the number of possible ways to iterate the equa- 
tions increase, creating a danger to run into a wrong so- 
lution. We finally explore other approximations to the 
exact solution that might be transposed to the full func- 
tional framework, namely a continued fraction approxi- 
mation and the expansion for large interaction, and we 
relate these approximations to the corresponding manip- 
ulations of the differential equation which produce them. 
These links arc crucial to prepare a generalization of the 
approach to the full functional framework. 
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Appendix A: Approximation for the Hartree term 

Due to the Hartree potential Vh = —ivG the set of dif- 
ferential equations (5) is nonlinear. In order to simplify 
this problem we first assume that Vh is Taylor expand- 
able in terms of the external potential tp: 



V H (3; [tp]) « -i / dM3+,4)G(4,4+; [tp]) 



<p=0 



i / d4d5v{3+A) 



«5G(4,4+;M) 



M5) 



<p(5) 



+ o(tp 2 ). (Al) 
The second step is to introduce G H defined through 
G°,(1,2) = G (1,2) 

d3G (l,3)V° (3)6^(3, 2), (A2) 



with V H {3) := -i J d4u(3+4)G(4,4 + ;fej]) 
ing Vh in Eq. (5) one obtains 



Insert- 



ed 



G(l, 2; [tp]) = G°,(l, 2) + / d3d5G H (l, 3) 
- i [ d4d5w(3 + ,4)- 



, ^8G{A,4+;[tp\) 



8tp(5) 



ip=0 



-5(3, 5)J tp(5)G(3,2;[tp]) 
-i J d3dAG H (l,3)v(3 + ,4) 



<K?(3,2;[y]) 
Stp(4) 



(A3) 



SG 

Since — — in the second term on the right-hand side of 

dip 

Eq. (A3) is a contraction of the two-particle correlation 
function, it yields the inverse dielectric function 



-i J d4v(3+,4) 
and one gets 



8G(4,A+;[tp}) 



8tp(5) 



ip=0 



+ 5(3,5) = e _1 (3, 5), 
(A4) 



G(l,2;[tp]) = G H (l,2) + J d3d5G° H (l,3) 

Xe _1 (3, 5)tp(5)G(3, 2; [<p]) + i J d3d4G° H (l,3) 

8G(3,2;[tp}) 
XV[6 ' 4) 8tp{4) ■ 



(A5) 



Now a rescaled perturbing potential can be introduced: 
tp(3) := / d5e- 1 (3,5)( /3 (5), (A6) 



, , , 8G 8G8<p . , , 
and, using the chain rule -j— = — — m the last term 

dip dip dip 

of the right-hand side of Eq. (A5), we get: 

G(l,2;fe>])=G&(l,2) + J d3d5G° H (l,3) 

xtp(3)G(3,2;[0\)+i J d3d5G%{l,3) 
8G(3,2;[tp}) 



xW(3 + ,5)- 



6<p(5) 



(A7) 



which is precisely Eq. (9). Here W = e 1 v is the screened 

Coulomb potential at vanishing (p. If one approximates 

8G 8G — ^ 

the functional derivative = —G — r^G « GG which 

dip dip 

comes from assuming the self-energy in the Dyson equa- 
tion G _1 = Gq" 1 — v H — £ — tp to be independent of tp, 
Eq. (A7) becomes: 

G(l,2;[tp])=G H (l,2) + / d3d5G° H (l,3) 



(A8) 



xtp(3)G(3,2;[<p}) + J d3d5G° H (l,3) 
xX GW (3,5;[tp})G(5,2;[tp]) 



with E G iv(3,5; [tp]) = £G(3,5; [tp])W(3 + ,5). For tp = 
Eq. (A8) becomes the Dyson equation for the one-particle 
Green's function in the GW approximation to the self- 
energy 35 . This confirms that the linearization of Vh is a 
reasonable starting point for further developments. 



Appendix B: Solving the DE 

Eq. (11) can be solved using standard textbook 
methods 36,37 . Here we choose a route that yields pre- 
cious information for our final aim of generalizing to the 
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full functional problem. A general ansatz for the struc- 
ture of y u {x) is: 



y u {x) = A(x) -l(x), 



(Bl) 



where the only restriction is that A and I are not zero. 
Substituting the ansatz in the DE (11) gives: 



dA(x) 

A(x)l(x) = y a + y a xA(x)2(x) - uy — — l(x) 

ax 



uy Q A(x 



dljx) 
dx 



(B2) 



The idea is now to solve two separate, simpler with re- 
spect to the initial one, DEs for A(x) and I{x). Putting 
together the left-hand side and the second and third 
terms of the right-hand side of Eq. (B2) one obtains: 

A{x)l{x) = y xA(x)l(x) - uy ^j^2{x). (B3) 



dx 



We can choose the solution 
A(x) = e 



(B4) 



which will then determine T(x) . One is now left with the 
equation for X(x) reading 



2/0 - uy Q A{x 



dl(x) 
dx 



0. 



(B5) 



Plugging in the expression for A(x) previously obtained 
and integrating on both sides one obtains: 



1 



dte^ 



T{x) 

u 

The integral on the right-hand side is: 

d 2 " + "™J = V2ue 2 " B o 

v/2u7T 



(B6) 



dt, 



uy ° dte 



x erf 



1 



yo 



/ r\ 2 



(B7) 



has 



where the change of variables t = ^-y== 

been made, and the lower limit of the last integral has 
been chosen to be zero, which requires to set a constant 
C(u,yo). Hence 



X(,) = y-e-o er/ 
+ C(u,y ). 



1 

yo 



l 



>2u 



(B8) 



The exact solution y u (x) = A(x) ■ T(x) is given in Eq. 

Flu - ~\, 

(14), where C{u,y ) = -\ — C(u, y )e 2 ^> . 



Appendix C: N-points continued fraction 
approximation 

We detail here how we have obtained the result of 
Eq. (53), or the order 0(d^.) of the N-points limited order 
DE. 

The starting point is Eq. (50), which is differentiated with 
respect to the external potential, yielding: 

5 2 G{l,2M) _ o C ^G(6,2; [y]) 

MWT^ fl(1 ' 6) " 



+ G° H (1J) 



8(p{7) 
8G(7,2; [0\) 



8<p(6) 



d3G° H (l,3)<p(3) 



8 2 G(3,2;m) 



^(6)^(7) 
i J d3d5W(3 + ,5)G° H (l,3) 
<J 3 G(3,2;[^]) 



Neglecting the term 
limit tp = yields: 

<5 2 G(1,2) 
^(7)^(6) = 



6<p(7)6<p(6)6(p(5)' 
^ 3 G(3,2;[^]) 

8G{&,2-[0. 



(CI) 

and taking the 



G° H (1,6) 
G° H (l,7) 



5<p(J) 
SG(7,2-M 



(p=0 



ip=0 



(C2) 



By substituting back Eq. (C2) into Eq. (50) we get: 
SG{1,2;[0\) 



ip=0 



G° H (1,6)G{6,2) 



+ i 



i J d3d5G° H (l,3)W(3 + ,5) 
8G(6, 2; [0\ 



G° H (3,6) 
G° ff (3,5) 



5tp(5) 
8G(5,2; [0\) 



<p=0 



8<p(6) 



tp=0 



(C3) 



The above equation can be recast in a compact way 

B xy = B° y + ^ l(xy)(qp)B qp , (C4) 



namely 



8G(l,2;[<p}) 



8<p(6) 



<p=0 



= 0^(1,6)0(6,2) 



* / d3d5d7G° H (l,3)W(3 + ,5)G° H (3,6) 
8G(7,2; [£ 



x 8(7,6)- 



8<p(5) 



ip=0 



d3d5d7G° H (l,3) 



x W(3 + ,5)G° H (3,b)8(7,6) 



SG(5,2-M 



8<p{7) 



tp=0 



(C5) 
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In the second term on the r.h.s. one can exchange, under 
the integral symbol, the indices 5 and 7, to obtain: 



5G(1,2; [0\) 



ip=0 



G° H (1,6)G(6,2) 

i J d3d5d7G° H {l,3)W(3 + ,5)G° H (3,6) 
5G(5,2;[0\) 



x 5(7,6)- 



ip=0 



x 5(7,6) 



d3d5d7G° H (l, 3)VF(3+, 5)G^(3, 5) 
5G(5,2;[v? 



ip=0 



(C6) 



Let's now define the following quantities: 



5G(1,2; [0\) 



Sq>(6) 

SG(5,2;[(p] 



6<p(7) 



ip=0 



ip=0 



:=.9(1,6), 



9(5,7), 



G° H {1,6)G(6,2) := 30(1,6), 

m(l,6;5,7) := i J d3G° H (l, 3)IT(3 + , 5)5(7, 6) 

[G° H (3,6) + G° H (3,5)]. (C7) 

Recasting Eq. (C6) with the new variables yields: 

#(1,6) =ff°(l,6) + I d5d7m(l,Q;5,7)g(5,7) (C8) 



we solve for g: 



d5d7 



m(l, 6; 5, 7) - 5(1, 5)5(7, 6) g(5,7) 

+ 5o(l,6) = 0. (C9) 

Defining 

m(l, 6; 5, 7) - 5(1, 5)5(7, 6)1 = m(l, 6; 5, 7), (CIO) 
inserting this expression into Eq. (C9): 

J dM7m{l, 6; 5, 7)g{h, 7) + g (l, 6) = 0, (Cll) 
and introducing the inverse of fh one gets: 

dld6d5d7m~ 1 (8, 9; 1, 6)m(l, 6; 5, 7)5(5, 7) = 
5m -1 (8, 9; l,6)flro(l, 6) (C12) 
5(8,9) = - / dld6m _1 (8, 9; 1,6)50(1, 6). (C13) 



Transforming back to the original variables gives 
8G(l,2;[0\) 



6<p(6) 



ip=0 



d9d8m _1 (l, 6; 9, 8)G^(9, 8)G(8, 2) 

(C14) 



and finally the Green's function reads: 
G{1,2) = G%{1,2) 

-i J d5d3d8d9G^(l,3)W(3 + ,5) 

x m- 1 (3,5;9,8)G^(9,8)G(8,2) (C15) 
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FIG. 1. Hedin's pentagon when W is kept fixed: one iterates 
only three equations, namely the ones for G, E, V rather than 
the full set. Note that keeping W fixed also implies fixing the 
polarizability P. 




FIG. 2. Comparison between the exact solution (red plain 
line, Eq. (15)) and the iterative solution for x — of 
(Eq. (20)). The blue crosses represent the I s * order expansion 
(Eq. (21)), while the green triangles and the black circles are 
respectively the 2 nd (Eq. (22)) and 3 rd order (Eq. (19)). All 
the three orders are close to the exact solution for small u val- 
ues, whereas when a given order of the series starts to diverge, 
the lower orders of the expansion reproduce the exact results 
better. For each curve C(u,yo) = —1, and we arbitrarily set 
yo = 1. 
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FIG. 3. Comparison between the exact solution (red plain 
line, Eq. (15)) and different flavours of the GW approxima- 
tion. In general the self-energy based approximations perform 
better than the iteration of the DE shown in Fig. 2. In the 
main panel the self-consistent GWo (black stars, Eq. (32)) is 
the best approximation to the exact result. Iterations start- 
ing from G — Go converge towards the self-consistent result 
(the 2 nd iteration is represented with light blue triangles, the 
3 rd with green circles and the 4 th with grey empty triangles). 
However, analyzing a larger u range (inset), one observes that 
odd iterations approach the exact u — oo limit, while the even 
ones don't seem to. It can be shown that they also do, how- 
ever in a very slow fashion and according to the following 
sequence t/i^L = {1/2, 1/3, 1/4, 1/5, 1/6, ■ ■ ■ }. 
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FIG. 4. In the main panel a comparison between the DE's 
exact solution (red plain line, Eq. (15)), GoW / oF' 1 - ) (blue 
squares, Eq. (42)), GWor' 1 ' (green empty circles, Eq. (43)) 
and sc-GWo (black stars, Eq. (32)) is shown. In this range 
of u, adding a vertex correction, no matter if within a self- 
consistent scheme or not, improves over the simpler self- 
consistent GWo- However, analyzing a wide it range (inset, 
semi- logarithmic plot), gives a different perspective: the first 
iteration of GqWoT^ 1 ' clearly exhibits the wrong u — > oo limit 
and the sc-GWo scheme becomes the closest approximation to 
the exact result. 
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FIG. 5. Comparison between the exact solution (red plain 
line, Eq. (15)) of the DE and the results obtained through the 
first three orders of the limited order differential equation (Eq. 
(49)). The notation C>(d™) indicates that derivatives of order 
> n have been neglected. As expected the result improves 
when more terms are included: the curve 0{di) (light blue 
line, Eq. (46)) is superimposed to the GoWo one (dark blue 
dots, Eq. (28)) and the curve 

0(4) (black circles, Eq. (48)) 
is close to the exact result in a small u range. 
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FIG. 6. Comparison between the exact solution (red plain 
line, Eq. (15)) and the large u expansion for the DE. The 
green stars and the black triangles are respectively 0(u _1//2 ) 
and 0(u~ 3 / 2 ) of the large u expansion (Eqs. (57-59)). We also 
report the GoWo result (blue dots, Eq. (28)) as an example 
of a small u expansion. Over a wide u range the large u 
expansions are very satisfactory. 
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FIG. 7. Comparison between the exact solution (plain red 
line, Eq. (15)), the G W result, the order 0{u~ 1/2 ) for the 
DE's large u expansion (green stars, Eq. (57) ), and the same 
order of the DE's large u expansion combined with the large 
£ expansion (black dots, Eq. (66)). We observe that the lat- 
ter approximation performs extremely well over the range of 
interaction examined, being even exact both in the large and 
small u limits. 



